home *** CD-ROM | disk | FTP | other *** search
/ Aminet 33 / Aminet 33 - October 1999.iso / Aminet / misc / math / TCalcStats2c.lha / TCalcStats2c / AREXX / Durbin.rexx < prev    next >
Encoding:
OS/2 REXX Batch file  |  1999-08-22  |  9.9 KB  |  519 lines
  1. /* Durbin Test*/
  2.  
  3. options results
  4. if ~show('P','TCALC') then do
  5.     address command 'run turbocalc:turbocalc'
  6.     address command 'waitforport TCALC'
  7.     loadflag=1
  8. end
  9. address 'TCALC'
  10. 'DEFPUBSCREEN()'
  11. /* Add-in Rexx Libraries needed for some routines */
  12. signal on syntax
  13. if ~show('l','rexxmathlib.library') then
  14.    call addlib('rexxmathlib.library',0,-30)  /* add to library list */
  15. if ~show('l','rexxreqtools.library') then
  16.    call addlib('rexxreqtools.library',0,-30)
  17. if ~show('l','rexxsupport.library') then
  18.    call addlib('rexxsupport.library',0,-30)
  19. signal off syntax
  20.  
  21. /* Start Main Routine */
  22. if loadflag=1 then 'Load()'
  23. 'ActivateWindow()'
  24. range=rtgetstring(,"Enter Cell Range for Input","Input Request") /*,,'rt_pubscrname="TCALC"')*/
  25. colon=pos(":",range)
  26. if colon=0 then do
  27.     'Message "Please select a range before executing this script"'
  28.     'DEFPUBSCREEN("Workbench")'
  29.     exit
  30. end
  31.  
  32. /* Find cell references and cell, column numbers */
  33. start_cell=substr(range,1,colon-1)
  34. end_cell=substr(range,colon+1)
  35. start_row=cellrow(start_cell)
  36. end_row=cellrow(end_cell)
  37. start_col=cellcol(start_cell)
  38. end_col=cellcol(end_cell)
  39. NRows=end_row-start_row+1
  40. NCols=end_col-start_col+1
  41.  
  42. /* Get cell reference for output range */
  43. out_cell=rtgetstring(,"Enter Cell Reference for Output","Input Request")/*,,'rt_pubscrname="TCALC"')*/
  44. if out_cell="" then exit
  45. if length(out_cell)<2 | datatype(left(out_cell,1),'n')=1 then do
  46.     'Message "Invalid cell reference"'
  47.     'DEFPUBSCREEN("Workbench")'
  48.     exit
  49. end
  50.  
  51. /* Suppress Screen Redraw to Speed Things Up */
  52. 'Refresh 0'
  53.  
  54. /* Open a small output window on tcalc screen*/
  55. fo=0
  56. CR='0a'x
  57. DisplayMsg="Calculating...Please Wait."||CR||"User input is disabled during calculation."||CR
  58. if open(6Info, 'con:100/0/450/80/Progress/SCREEN TCALC', w) then do
  59.      call writeln(6Info, DisplayMsg)
  60.     fo=1
  61. end
  62. else do
  63.     'Message "TCALC Screen not available for Progress messages"'
  64. end
  65. CALL DELAY(150)
  66.  
  67. /* Get cell references for top cell in each column */
  68. 'SelectCell' start_cell
  69. do col=start_col to end_col
  70.     'GetCursorPos'
  71.     top_cell.col=result
  72.     'Column' 1
  73. end
  74.  
  75. /* Get labels for later use on output */
  76. 'SelectCell' start_cell
  77. 'GetValue'
  78. testlabel=result
  79. testlabel=strip(testlabel)
  80. if datatype(testlabel,'n')=1 then do
  81.     labelflag=0
  82.     do x=1 to NCols
  83.     title.x="Column "||x
  84.     end
  85. end
  86. else do
  87.     labelflag=1
  88.     NRows=NRows-1
  89.     do x=1 to NCols
  90.     'GetValue'
  91.     str=result
  92.     title.x=translate(strip(str),"_"," ")
  93.     'Column 1'
  94.     end
  95. end
  96. if fo then call writech(6Info,"Progress...10 ")
  97.  
  98. /* Get data from cell range */
  99. col=start_col
  100. lav=0
  101. tot=0
  102. count.=0
  103. R.=0
  104. do x=1 to NCols
  105.     'SelectCell' top_cell.col
  106.     if labelflag=1 then 'CursorDown 1'
  107.     do y=1 to NRows
  108.         'GetValue'
  109.         valtest=result
  110.         if datatype(valtest)='NUM' then do
  111.             'GetValue'
  112.             val=result
  113.             data.x.y=val
  114.             tot=tot+val
  115.             R.x=tot
  116.             count.x=1+count.x
  117.         end
  118.         'CursorDown 1'
  119.     end
  120.     col=col+1
  121.     tot=0
  122.     lav=0
  123.     val=0
  124. end
  125.  
  126. if fo then call writech(6Info,"20 ")
  127.  
  128. /* Are there any ties? */
  129. temp=0
  130. Ties=0
  131. DO y=1 to NRows
  132.     Do x=1 to NCols
  133.         If datatype(data.x.y)='NUM' then temp=data.x.y
  134.         ELSE iterate
  135.         Do z=1 to NCols
  136.             BZ=0
  137.             IF datatype(data.z.y)='NUM' then do
  138.                 IF temp=data.z.y then BZ=BZ+1
  139.                 End
  140.         IF BZ>1 then Ties=1
  141.         End
  142.     End
  143. End
  144. if fo then call writech(6Info,"40 ")
  145. say "Ties="Ties
  146. /* Number of values in each row and column? */
  147. r=count.1
  148. k=0
  149. y=1
  150. Do x=1 to NCols
  151.         if datatype(data.x.y)='NUM' then k=k+1
  152. End
  153. /* Calculate T1 and T2*/
  154. T1=0
  155. T2=0
  156. Ta=0
  157. Tb=0
  158. Tc=0
  159. Td=0
  160. IF Ties=0 then Do
  161.     Ta=(12*(NCols-1))/(r*NCols*(k-1)*(k+1))
  162.     Do x=1 to NCols
  163.         Tb=Tb+(R.x-(r*(k+1))/2)**2
  164.     End
  165.     T1=Ta*Tb
  166.     Tc=T1/(NCols-1)
  167.     Td=(NRows*(k-1)-T1)/(NRows*k-NRows-NCols+1)
  168.     T2=Tc/Td
  169. End
  170. Else Do
  171.     A=0
  172.     C=0
  173.     Te=0
  174.     Do  x=1 to NCols
  175.         Do y=1 to NRows
  176.             if datatype(data.x.y)='NUM' then A=A+(data.x.y)**2
  177.         End
  178.     End
  179.     C=(NRows*k*(k+1)**2)/4
  180.     Do x=1 to NCols
  181.         Te=Te+(R.x-(r*(k+1))/2)**2
  182.     End
  183.     Ta=(NCols-1)*Te
  184.     Tb=A-C
  185.     T1=Ta/Tb
  186.     Tc=T1/(NCols-1)
  187.     Td=(NRows*(k-1)-T1)/(NRows*k-NRows-NCols+1)
  188.     T2=Tc/Td
  189. End
  190. if fo then call writech(6Info,"60 ")
  191.  
  192. /* Calculate Probability */
  193.  
  194. AC=.0001
  195. EC=.005
  196. EC2=EC+EC
  197. DF1=NCols-1
  198. DF2=NRows*k-NRows-NCols+1
  199. P=PROBABILITY(T2,DF1,DF2)
  200. P=1-P
  201. P2=P*2
  202. if fo then call writech(6Info,"80 ")
  203. /* P=.95 OR .99 */
  204. FCRIT1=FCRITICAL(.95,DF1,DF2)
  205. FCRIT2=FCRITICAL(.99,DF1,DF2)
  206. FCRIT3=FCRITICAL(.975,DF1,DF2)
  207. FCRIT4=FCRITICAL(.995,DF1,DF2)
  208.  
  209. if fo then do
  210.     call writeln(6Info,"100 ")
  211.     call writeln(6Info,"Writing output to window...")
  212. end
  213.  
  214. /* Output */
  215. 'SelectCell' out_cell
  216. 'ColumnWidth 18'
  217. 'Put "Durbin Test - Balanced Incomplete Block Design"'
  218. 'CursorDown 1'
  219. IF Ties=0 then 'Put "(There were no ties)"'
  220. ELSE 'Put "(Ties were detected)"'
  221. 'CursorDown 1'
  222. 'Put "Durbin T:"'
  223. 'CursorDown 1'
  224. 'Put "df(1):"'
  225. 'CursorDown 1'
  226. 'Put "df(2):"'
  227. 'CursorDown 1'
  228. 'Put "P(F<=f) One-tail:"'
  229. 'CursorDown 1'
  230. 'Put "F-Critical(95%):"'
  231. 'CursorDown 1'
  232. 'Put "F-Critical(99%):"'
  233. /*'CursorDown 1'
  234. 'Put "P(F<=f) Two-tail:"'
  235. 'CursorDown 1'
  236. 'Put "F-Critical(95%):"'
  237. 'CursorDown 1'
  238. 'Put "F-Critical(99%):"'*/
  239. 'SelectCell' out_cell
  240. 'CursorDown 2'
  241. 'CursorRight 1'
  242. 'Put' format(T2,,4)
  243. 'CursorDown 1'
  244. 'Put' DF1
  245. 'CursorDown 1'
  246. 'Put' DF2
  247. 'CursorDown 1'
  248. 'Put' format(P,,6)
  249. 'CursorDown 1'
  250. 'Put' format(FCRIT1,,4)
  251. 'CursorDown 1'
  252. 'Put' format(FCRIT2,,4)
  253. /*'CursorDown 1'
  254. 'Put' format(P2,,6)
  255. 'CursorDown 1'
  256. 'Put' format(FCRIT3,,4)
  257. 'CursorDown 1'
  258. 'Put' format(FCRIT4,,4)*/
  259. 'Refresh 1'
  260. 'Refresh 2'
  261. /*'Message' "Finished"*/
  262. /*indicate the main script is finished*/
  263. DisplayMsg="Cleaning up ...."||CR||"Exiting"
  264. result=writeln(6Info, DisplayMsg)
  265. if result~=0 then do
  266.     /*Wait 3 seconds*/
  267.     CALL DELAY(150)
  268.     /* close window*/
  269.     result=close(6Info)
  270. end
  271. 'DEFPUBSCREEN("Workbench")'
  272. exit
  273.  
  274. /* Procedures */
  275.  
  276. cellrow: procedure
  277. do
  278.     parse arg cell
  279.     do charpos=2 to length(cell)
  280.     if datatype(substr(cell,charpos,1),n) then return substr(cell,charpos)
  281.     end
  282.     return 0
  283. end
  284. Return
  285.  
  286. cellcol: procedure
  287. do
  288.     parse arg cell
  289.     labels="ABCDEFGHIJKLMNOPQRSTUVWXYZ"
  290.     cell=upper(cell)
  291.     len=length(cell)
  292.     val=0
  293. do charpos=1 to len
  294.     if datatype(substr(cell,charpos,1),n) then
  295.     do cell=reverse(substr(cell,1,charpos-1))
  296.     do x=1 to length(cell)
  297.     val=(26**(x-1))*pos(substr(cell,x,1),labels)+val
  298.     end
  299.     return val
  300.     end
  301.     end
  302.     return 0
  303. end
  304. Return
  305.  
  306. syntax:
  307.      if arg(1)='FAIL' then do
  308.         'Message "Library is unavailable."'
  309.         'DEFPUBSCREEN("Workbench")'
  310.         exit
  311.         end
  312.     'Message "Unknown Syntax Error"'
  313.     'DEFPUBSCREEN("Workbench")'
  314.     exit
  315.  
  316. Format:  procedure
  317.  
  318.      arg number, before, after
  319.      CallLine = SIGL
  320.      if ~datatype(CallLine, 'N') then CallLine = '??'
  321.  
  322.      /* Make sure we have a number as first (required) argument    */
  323.      if ~datatype(number, 'N') then do
  324.         if number = '' then
  325.            fc = 17     /* Wrong number of arguments           */
  326.         else
  327.            fc = 47     /* Arithmetic conversion error             */
  328.         signal FormatSyntaxError
  329.      end
  330.      num = number + 0
  331.      if before = '' & after = '' then
  332.         return num
  333.      else do
  334.         parse var num integer '.' fraction
  335.         if before = '' then before = length(integer)
  336.         if after = '' then after = length(fraction)
  337.         if ~datatype(before, N) | ~datatype(after, N) then
  338.            do fc = 18
  339.            signal FormatSyntaxError
  340.        end
  341.         if before < length(integer) then do
  342.            fc = 18
  343.            signal FormatSyntaxError
  344.         end
  345.         if after ~= length(fraction) then do
  346.            fraction = trunc(('.'fraction'0') + ('.'copies('0', after)'5'), after)
  347.         if integer<1&integer>-1 then integer=integer
  348.            else integer = integer + (fraction % 1)
  349.            fraction = substr(fraction, 3)
  350.         end
  351.         if fraction >= 0 then
  352.            return right(integer, before)'.'fraction
  353.         else
  354.            return right(integer, before)
  355.      end
  356.  
  357.  FormatSyntaxError:
  358.         if show('F', STDERR) then
  359.            call writeln(STDERR, '+++ Error' fc 'in line' CallLine':' errortext(fc))
  360.         else
  361.            mess='+++ Error' fc 'in line' CallLine':' errortext(fc)
  362.         'Message' mess
  363.         parse source Func .
  364.         if Func = 'FUNCTION' then do
  365.         'DEFPUBSCREEN("Workbench")'
  366.            exit "Err"
  367.         end
  368.         else do
  369.         'DEFPUBSCREEN("Workbench")'
  370.            exit 10
  371.         end
  372.  
  373. FCRITICAL: PROCEDURE EXPOSE AC EC EC2
  374.  
  375.     PARSE ARG PO,K1,K2
  376.     /* fIND NORMAL Z CORRESPONDING TO P */
  377.     P1=PO
  378.     Z=NORMALPP(PO)
  379.     IF K2<3 THEN DO
  380.         W=K2**.75
  381.         W=Z/W
  382.         Z=Z*(1.1581-W*(.2296+W*(.0042+.0027*W)))
  383.     END
  384.     /* FIND INITIAL APPROX. TO F */
  385.     A1=2/(9*K1)
  386.     A2=2/(9*K2)
  387.     W=Z*Z
  388.     W1=1+A2*(A2-W-2)
  389.     W2=A1+A2-A1*A2-1
  390.     W3=1+A1*(A1-W-2)
  391.     SN=SIGN(W2*W2-W1*W3)
  392.     IF SN=-1 THEN W3=-1*SQRT(ABS(W2*W2-W1*W3))
  393.     ELSE W3=SQRT(W2*W2-W1*W3)
  394.     FO=(W3-W2)/W1
  395.     FO=FO**3
  396.     /* MODIFIED NEWTON ITERATION TO IMPROVE VALUE OF F */
  397.     DO FOREVER
  398.         FCRIT=FO
  399.         PO=PROBABILITY(FCRIT,K1,K2)
  400.         F1=PO-P1
  401.         FCRIT=FO+EC
  402.         PO=PROBABILITY(FCRIT,K1,K2)
  403.         F2=PO
  404.         FCRIT=FO-EC
  405.         PO=PROBABILITY(FCRIT,K1,K2)
  406.         F2=F2-PO
  407.         F2=F2/EC2
  408.         F1=FO-F1/F2
  409.         IF ABS(F1-FO)>AC THEN FO=F1
  410.         ELSE LEAVE
  411.     END
  412.     FCRIT=F1
  413. RETURN FCRIT
  414.  
  415. NORMALPP: PROCEDURE
  416.  
  417.     ARG P0
  418.     A1=2.515517
  419.     A2=.802853
  420.     A3=.010328
  421.     A4=1.432788
  422.     A5=.189269
  423.     A6=.001308
  424.     Q=.5-ABS(P0-.5)
  425.     SN=SIGN(-2*LN(Q))
  426.     IF SN=-1 THEN W=-1*SQRT(ABS(-2*LN(Q))
  427.     ELSE W=SQRT(-2*LN(Q))
  428.     W1=A1+W*(A2+A3*W)
  429.     W2=1+W*(A4+W*(A5+A6*W))
  430.     ZZ=W-W1/W2
  431.     SELECT
  432.         WHEN (P0-.5)<0 THEN TT=-1
  433.         WHEN (P0-.5)=0 THEN TT=0
  434.         otherwise TT=1
  435.     END
  436.     ZZ=ZZ*TT
  437. RETURN ZZ
  438.  
  439. PROBABILITY: PROCEDURE EXPOSE AC EC EC2
  440.  
  441.     PARSE ARG F,K1,K2
  442.     IF K1=1 THEN AC=.00001
  443.     H2=K2/2
  444.     H1=K1/2
  445.     H3=H1+H2
  446.     L1=0
  447.     XX=H1
  448.     LG=LOGGAMMA(XX)
  449.     L1=L1+LG
  450.     XX=H2
  451.     LG=LOGGAMMA(XX)
  452.     L1=L1+LG
  453.     XX=H3
  454.     LG=LOGGAMMA(XX)
  455.     L1=L1-LG
  456.     W1=K2/(K2+K1*F)
  457.     XX=0
  458.     IF H2<(H3*W1) THEN DO
  459.         W2=W1
  460.         W1=1-W1
  461.         W3=H2
  462.         H2=H1
  463.         H1=W3
  464.     END
  465.     ELSE DO
  466.         W2=1-W1
  467.         XX=1
  468.     END
  469.     T1=1
  470.     W3=1
  471.     P=1
  472.     I=INT(H1+W2*H3)
  473.     W4=W1/W2
  474.     /*LABELB:*/
  475.     T2=H1-W3
  476.     IF I=0 THEN W4=W1
  477.     /*LABELC:*/
  478.     DO FOREVER
  479.         T1=T1*T2*W4/(H2+W3)
  480.         P=P+T1
  481.         T2=ABS(T1)
  482.         IF (T2<=AC) & (T2<=AC*P) THEN LEAVE /*SIGNAL 'LABELD'*/
  483.         W3=W3+1
  484.         I=I-1
  485.         IF I>=0 THEN DO /*SIGNAL 'LABELB'*/
  486.             T2=H1-W3
  487.             IF I=0 THEN W4=W1
  488.             ITERATE
  489.         END
  490.         T2=H3
  491.         H3=H3+1
  492.     /*SIGNAL 'LABELC'*/
  493.     END
  494.     /*LABELD:*/
  495.     W1=EXP(H2*LN(W1)+(H1-1)*LN(W2)-L1)
  496.     P=P*W1/H2
  497.     IF XX=0 THEN PO=P
  498.     ELSE PO=1-P
  499.     AC=.001
  500. RETURN PO
  501.  
  502. LOGGAMMA: PROCEDURE
  503.  
  504.     ARG XA
  505.     C1=76.18009173
  506.     C2=-86.50532033
  507.     C3=24.01409822
  508.     C4=-1.231739516
  509.     C5=.001208580
  510.     C6=-.000005364
  511.     C7=2.506628275
  512.     X1=XA-1
  513.     W=X1+5.5
  514.     W=(X1+.5)*LN(W)-W
  515.     S=1+C1/(X1+1)+C2/(X1+2)+C3/(X1+3)
  516.     S=S+C4/(X1+4)+C5/(X1+5)+C6/(X1+6)
  517.     L=W+LN(C7*S)
  518. RETURN L
  519.